TODO:
Calculate the upwind area aka footprint (\(X_f\), m) representative of these measurements that accounts for a percentage of the flux. Where \(U\) is the assumed constant wind speed, \(u^{*}\) is the friction velocity, \(d\) is the displacement height due to vegetation (taken as 2/3 the vegetation height), \(z\) is the instrument height, \(k\) is the von karmans constant (usually taken as 0.4), \(f\) is the fraction of the area that explains the flux [@Schuepp1990]:
\[ X_f = - \frac{1}{ln(f)}\frac{U}{u^*}\frac{z-d}{k} \] Friction velocity (\(u_*\)) can be modelled as in wind_profile.Rmd or is also measured at the Eddy Covariance Sites.
Constants:
# need to measure this in the field still (2022-06-01)
temp_offset <- 125
k <- 0.4 # von karmans constant (unitless)
d0 <- 1.965869 # m estimated from jan 16 2022 at 4:00 pm
f <- 0.8 # 80%
Load Data:
low_wind <- readRDS('data/eddy_cov_15min_combined.rds') |> select(TIMESTAMP, wnd_dir_compass, u_star, low_ec_rslt_wnd_spd = result_wind_spd)# |>
# pivot_longer(u_star:low_ec_rslt_wnd_spd)
#
# ggplot(low_wind, aes(TIMESTAMP, value, colour = name)) + geom_line()
#
# plotly::ggplotly()
mid_hi_wind <- readRDS('../../analysis/treefort-main/data/met/met_main.rds') |> select(TIMESTAMP, USWindSpeed_S_WVT, top_ec_rlst_wnd_spd, top_ec_u_star)
wind_wide <- left_join(low_wind, mid_hi_wind, by = 'TIMESTAMP')
wind_dir <- wind_wide |>
mutate(wind_dir = (wnd_dir_compass + temp_offset) %% 365) |>
select(TIMESTAMP, wind_dir)
low <- wind_wide |>
select(TIMESTAMP, u_star, wind_speed = low_ec_rslt_wnd_spd) |>
mutate(inst_height = 3) |>
mutate(
med = zoo::rollapply(wind_speed, width = 1000, FUN = median, na.rm = T, align = 'center', partial = T),
sd = zoo::rollapply(wind_speed, width = 1000, FUN = sd, na.rm = T, align = 'center', partial = T),
stdep = (wind_speed - med) / sd,
flag = ifelse(wind_speed > 4, T, F)
) |>
select(TIMESTAMP, u_star, wind_speed, inst_height, flag)
ggplot(low, aes(TIMESTAMP, wind_speed, colour = flag)) + geom_point()
## Warning: Removed 1372 rows containing missing values (geom_point).
plotly::ggplotly()
hi <- wind_wide |>
select(TIMESTAMP, u_star = top_ec_u_star, wind_speed = top_ec_rlst_wnd_spd) |>
mutate(
inst_height = 13.5,
mean = zoo::rollapply(wind_speed, width = 48, FUN = mean, na.rm = T, align = 'center', partial = T),
sd = zoo::rollapply(wind_speed, width = 48, FUN = sd, na.rm = T, align = 'center', partial = T),
stdep = (wind_speed - mean) / sd,
flag = ifelse(wind_speed > 20, T, F)
) |>
select(TIMESTAMP, u_star, wind_speed, inst_height, flag)
ggplot(hi, aes(TIMESTAMP, wind_speed, colour = flag)) + geom_point()
## Warning: Removed 8863 rows containing missing values (geom_point).
plotly::ggplotly()
all <- rbind(low, hi) |>
left_join(wind_dir) |>
filter(flag == F)
## Joining, by = "TIMESTAMP"
all <- weatherdash::cat_wind(all, 'wind_dir', 'wind_dir_bin')
ggplot(all, aes(TIMESTAMP, wind_speed, colour = flag)) +
geom_point() +
facet_wrap(~inst_height)
plotly::ggplotly()
#|> filter(as.Date(TIMESTAMP) == as.Date('2022-01-16'))
Calculate \(X_f\):
x_f <- all |>
mutate(x_f = -(1/log(f))*(wind_speed/u_star)*((inst_height-d0)/k)) |>
filter(x_f < 5000)
Plot \(X_f\) for 3m EC with wind speed on all days: NOTE EC WAS RAISED Feb 17, 2022, this could be why relationship changes after FEB plot.
ggplot(x_f |> filter(inst_height == 3), aes(wind_speed, x_f, colour = wind_dir_bin)) +
geom_point() +
facet_wrap(~inst_height+strftime(as.Date(TIMESTAMP), "%m"))
Plot \(X_f\) for 13.5 m EC with wind speed on all days:
Maybe larger footprint during winter months when snow is in the trees?
ggplot(x_f |> filter(inst_height == 13.5), aes(wind_speed, x_f, colour = wind_dir_bin)) +
geom_point() +
facet_wrap(~inst_height+strftime(as.Date(TIMESTAMP), "%m"), scales = 'free')
plotly::ggplotly()
Output \(X_f\) table
x_f |>
mutate(yr_mon = strftime(as.Date(TIMESTAMP), "%y-%m")) |>
group_by(inst_height, yr_mon) |>
summarise(x_f_mean = mean(x_f, na.rm = T),
u_star_mean = mean(u_star, na.rm = T),
wind_spd_mean = mean(wind_speed, na.rm = T),
wind_spd_pk_mean = max(wind_speed, na.rm = T),
wind_dir_mean = mean(wind_dir, na.rm = T))
## `summarise()` has grouped output by 'inst_height'. You can override using the
## `.groups` argument.
## # A tibble: 14 × 7
## # Groups: inst_height [2]
## inst_height yr_mon x_f_mean u_star_mean wind_spd_mean wind_spd_pk_mean
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 3 21-11 35.4 0.226 0.613 2.88
## 2 3 21-12 42.8 0.186 0.524 2.85
## 3 3 22-01 55.6 0.186 0.796 2.82
## 4 3 22-02 52.2 0.167 0.629 2.62
## 5 3 22-03 46.4 0.162 0.557 3.74
## 6 3 22-04 53.4 0.141 0.510 1.72
## 7 3 22-05 51.3 0.189 0.588 2.23
## 8 13.5 21-11 1165. 0.500 4.15 10.1
## 9 13.5 21-12 1118. 0.450 3.53 9.17
## 10 13.5 22-01 1097. 0.524 4.02 11.1
## 11 13.5 22-02 1061. 0.469 3.39 9.12
## 12 13.5 22-03 957. 0.443 2.96 10.4
## 13 13.5 22-04 930. 0.420 2.75 9.56
## 14 13.5 22-05 1055. 0.442 3.21 8.71
## # … with 1 more variable: wind_dir_mean <dbl>
Wind Rose for all months for each instrument height:
# for (i in c(3,13.5)) {
# for (m in 1:12) {
# fltr <- all |>
# mutate(month = as.numeric(strftime(as.Date(TIMESTAMP), "%m"))) |>
# filter(inst_height == i,
# month == m)
#
# if(nrow(fltr) != 0){
# p <- weatherdash::wind_rose(fltr, 'TIMESTAMP', 'wind_speed', 'wind_dir', plot_title = paste('Height', i, "m", 'for', m))
# file1 <- paste0('figs/wind_rose/wind_rose_height_', i, '_', m,'.html')
# file2 <- paste0('figs/wind_rose/wind_rose_height_', i, '_', m,'.png')
#
# htmlwidgets::saveWidget(plotly::as_widget(p), file1)
# plotly::save_image(p, file2)
# }
#
# }
#
# }
Look at one of the windiest days without snow on instruments:
ggplot(x_f |> filter(as.Date(TIMESTAMP) == as.Date('2022-01-16')) , aes(wind_speed, x_f)) +
geom_point() +
facet_wrap(~inst_height, scales = 'free')
Wind Rose on this Day:
low_wind <- readRDS('data/eddy_cov_15min_combined.rds') |>
filter(as.Date(TIMESTAMP) == as.Date('2022-01-16')) |>
select(TIMESTAMP,
wnd_dir_compass,
u_star,
low_ec_rslt_wnd_spd = result_wind_spd) |>
mutate(
wind_dir = (wnd_dir_compass + temp_offset) %% 365,
mean = zoo::rollapply(low_ec_rslt_wnd_spd, width = 48, FUN = mean, na.rm = T, align = 'center', partial = T),
sd = zoo::rollapply(low_ec_rslt_wnd_spd, width = 48, FUN = sd, na.rm = T, align = 'center', partial = T),
stdep = (low_ec_rslt_wnd_spd - mean) / sd,
flag = ifelse(low_ec_rslt_wnd_spd > 4, T, F)
)
ggplot(low_wind, aes(TIMESTAMP, low_ec_rslt_wnd_spd, colour = flag)) +
geom_point()
weatherdash::wind_rose(low_wind, 'TIMESTAMP', 'low_ec_rslt_wnd_spd', 'wind_dir', plot_title = 'test')